, •• >c 

NASA Technical Memorandum 87799 


A Method of Calculating the 
Total Flow From a Given 
Sea Surface Topography 


Desiraju B. Rao, Stephen D. Steenrod, 
and Braulio V. Sanchez 


APRIL 1987 




NASA Technical Memorandum 87799 


A Method of Calculating the 
Total Flow From a Given 
Sea Surface Topography 


Desiraju B. Rao 

National Meteorological Center 

Washington, D.C. 

Stephen D. Steenrod 
Applied Research Corporation 
handover, Maryland 

Braulio V. Sanchez 
Goddard Space Flight Center 
Greenbelt, Maryland 


NASA 

National Aeronautics 
and Space Administration 

Scientific and Technical 
Information Branch 


1987 


TABLE OF CONTENTS 


INTRODUCTION 1 

THEORETICAL BASIS 3 

RESULTS FOR AN IDEAL CASE 7 

CONCLUSION 15 

REFERENCES 19 


PRECEDMQ page blank hot hlmed 


A Method of Calculating the Total Flow From a Given Sea Surface Topography 


Desiraju B. Rao, Stephen D. Steenrod, and Braulio V. Sanchez 

INTRODUCTION 

Altimetric measurements from a satellite yield the sea surface topography on a global basis. Such a sea 
surface topography represents the sum of (i) the marine geoid, which is the shape the ocean surface assumes 
when at rest, (ii) the departure of the ocean surface from the geoid due to the near surface ocean circula- 
tions, and (iii) a variety of errors related to the determination of satellite orbit, tides, etc. See Wunsch and 
Gaposchkin (1980) and report of the TOPEX Science Working Group (1981) for a review of the errors 
associated with satellite altimetric measurements. If the errors due to source (iii) are eliminated and if one 
knows the shape of the geoid, in principle, we have the residual sea surface topography associated with 
dynamics of the ocean circulation. Since the surface elevation is a simple scalar variable and reflects ocean- 
ic processes occurring at depth, satellite altimetry offers the potential to study global ocean circulation if 
it can be measured with the necessary accuracy - particularly when treated together with other available 
hydrographic data. See, for example, the analysis of the Seasat altimetric residuals by Tai and Wunsch 
(1983) and Tai (1983). 

The primary advantage of being able to measure the sea surface topography is the fact that one can 
directly compute the surface pressure gradients and the associated geostrophic currents. These currents would 
then be the reference velocities for the integration of the observed density field via the “thermal wind” 
equation to provide absolute currents at depth. Hence the altimetric measurements of the sea surface eliminate 
the necessity of invoking the existence of an arbitrary level of no motion which can only provide the relative 
currents. However, if certain assumptions are invoked and a conceptual model of ocean circulation dynam- 
ics is adapted, it is possible to infer more than the “near surface geostrophic currents” from the measured 
sea surface topography. Admittedly the ocean circulation is a complex result of forcing by the wind stresses 
and the atmospheric heating and cooling. The circulation consists of both barotropic and baroclinic compo- 
nents and exhibits spatial and temporal variabilities on a wide range of scales. Nevertheless, in this prelimi- 
nary attempt, a working assumption is made that the circulation dynamics in the upper layers of the ocean 
are dominated by wind driven forcing much like in the Stommel (1948) study. When the upper layer of 
an ocean is spun up to a steady state under the influence of a wind stress and a linear bottom friction, the 
circulation dynamics represent a balance between the Coriolis, pressure gradient, wind stress, and bottom 
friction forces. Hence the observed surface topography represents a more complicated balance of forces 
than just geostrophy even in this simple conceptual dynamical framework. 

We present a method here that, given a data field of sea surface topography, allows the recovery through 
an objective analysis procedure of the total flow field, which may be viewed as the sum of a geostrophic 
and an ageostrophic component if one so desires. However, since geostrophic considerations become progres- 
sively weak as one approaches the equator and the flows become more and more ageostrophic, we will 
simply refer to the flow derived from our calculations as the actual or total flow. Since no geostrophic con- 
siderations are imposed in the analytical development, the analysis procedure outlined here is not affected 
by the vanishing of the Coriolis parameter at the equator in calculating the currents from the surface pres- 
sure field. 

The results presented deal with the steady circulation in an ideal rectangular basin on a beta plane. The 
theory is, however, also valid to study the transient circulation so long as the circulation is assumed to 
be non-divergent. 

The theory is based on developing a spectrum of characteristic functions for the stream function field 
and the surface height field. Given a wind field these functions are then used to represent the forced solution 
for the velocity and surface pressure fields in which the expansion coefficients for both fields are required 
to be the same. Such a calculation serves a two-fold purpose. The first one is that the forced solution indi- 
cates which of the spectral components are most energetic so that these may be used in the objective analy- 
sis of the height field data. The second purpose is that when an objective analysis is performed on a limited 
set of surface height field data points extracted from the theoretical solution, this solution would serve the 
function of the “true” solution to validate the objective analysis procedure. The reason for using a limited 
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set of available data in dealing with real altimetric signals is simply due to the errors in our knowledge 
of the geoid over various parts of the global ocean and the consequent contamination of the altimetric retrievals 
over these regions of the global ocean. Elimination of data in regions of serious uncertainty in the altimetric 
retrievals and the use of functions characteristic of a given basin should result in extracting the oceano- 
graphic signal from the measurements in a most optimal way. In the objective analysis the most energetic 
functions, as indicated by the theoretical solution, would be used and the coefficients of expansion would 
be determined by linear least squares procedure. Once the expansion coefficients are thus obtained, the 
stream function associated with the total flow can be synthesized using the stream function characteristic 
modes. An additional advantage of the proposed technique is that the analytical basis of the theory allows 
inferences to be made on the direction of the rotational part of the wind stress. At least in the simple case 
presented here, the direction of the wind stress is unambiguously determined. The entire methodology 
described here is easily treated numerically and can be extended to deal with real world basins as will be 
discussed later. 
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THEORETICAL BASIS 


In solving the problem of wind driven ocean circulation 
in an enclosed basin using “spectral” expansion proce- 
dures, it is possible to consider the full gravitational and rota- 
tional modes of the basin as the relevant functions. The 
powerful expansion technique of Eckart may then be invoked 
to solve the inhomogeneous problem as, for example, done 
by Reid (1958) and Rao (1974). Although the determination 
of gravitational normal modes of a closed basin is not an easy 
matter even for simple shapes, it is nevertheless feasible (Rao 
1966). Moreover, even for natural basins numerical methods 
have now been developed to compute the gravitational modes 
reasonably accurately (Platzman 1975, 1977, Rao and 
Schwab 1976). Rotational modes, however, are very com- 
plex even in simple basins and for real basins they are ex- 
tremely sensitive to the resolution of the shape of the basins 
and the bottom topographic variations as found in the above 
mentioned studies of Platzman, and Rao and Schwab. Since 
the large scale ocean circulation is primarily dominated by 
rotational modes and the determination process of these 
modes is not a “robust” one, it is felt that it would be best 
to avoid a representation of the solution in terms of normal 
modes. Instead, the solutions are developed in terms of a 
set of characteristic functions that are more easily determina- 
ble and are less sensitive to the resolution of the basin ge- 
ometry and bottom topography. Further, this procedure 
avoids the dynamic “prejudice” such as absence of friction, 
sensitivity to spatial resolutions, etc., implicit in the actual 
normal modes. Indeed this is the procedure adapted by Rao 
and Schwab (1981) in the analysis of lake circulations and 
by Sanchez, Rao, and Wolfson (1985) in the analysis of tides 
in a lake. 

The analysis presented below depends primarily on the 
basic assumptions that the ocean circulation is non-divergent 
and that the wind stress acting on it is also non-divergent. 
For large scale ocean circulation these assumptions can be 
justified using the conventional scale considerations. An ad- 
ditional assumption introduced is that the normal component 
of the wind stress is zero on the north-south boundaries for 
reasons to be explained later. In this study we consider a 
linear wind driven ocean circulation model with a linear bot- 
tom friction. The analysis can, however, be extended to take 
into account nonlinear dynamics if necessary. The govern- 
ing equations are: 


du 

dt 


fv = 


8 S + 



Xu 


dv 

at 


+ 


fu = 


ar , II 

dy pH 


Xv 


( 2 . 1 ) 


In the above equations, u and v are the horizontal veloc- 
ities in the X (east-west) and Y (north-south) directions, f 
is the Coriolis parameter = f 0 + /Sy. f represents the sur- 
face height fluctuation (namely, the altimetric residual of the 
sea surface), g the gravitational acceleration, r x , r y the wind 
stress, and X the linear friction coefficient, q is the water 
density which is taken to be = 1 cgs and omitted from here 
on. The depth H in equation (2.1) is the so called depth of 
wind driven circulation. This depth will be taken as a con- 
stant even though such a restriction is not necessary for the 
analysis. If it is considered that the wind driven circulation 
is felt through the entire depth of the ocean and the depth 
H is a variable, the analysis can easily accommodate this fol- 
lowing the line of Rao and Schwab’s (1981) study. 

In view of the fact that H is a constant and the flow is 
non-divergent, the velocity field can be expressed in terms 
of a stream function: 


u = — — — v — — — 
dy ’ dx 


(2.3) 


Forming the vorticity equation from (2.1), we obtain 

(|+X) , + / 3v=i(^-^) (2.4) 

where the vertical component of vorticity is given by 
V = (2.5) 


If the divergence operator is applied to equation (2. 1) the 
result, in view of equation (2.2) and the assumption that the 
large scale wind stress is also non-divergent, is: 

v.f ViA = gV 2 ^ (2.6) 

The above equation is easily recognized as the “linear 
balance equation” extensively used by the meteorological 
community in numerical weather prediction. 

Even though equation (2.6) appears to be a direct conse- 
quence of the geostrophic relations 

f Vtf = gVf , 


it should be noted that exact geostrophy is not permitted for 
two dimensional nondivergent flows on a rotating earth with 
a variable Coriolis parameter. This can be seen by taking 
the curl of the above relations which leads to the result that 


9fB± = n 

dy dx 


du dv 
dx dy 


0 


( 2 . 2 ) 


Since df/dy =£ 0, then v = d\j//dx must be = 0. The con- 
tinuity equation then reduces to du /dx = 0. Hence u itself 
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is zero since it has to vanish on the eastern and western 
boundaries of the basin leading to the conclusion that non- 
divergent geostrophic flow is not permitted on a beta-plane 
or a spherical earth. On the other hand if the bottom topog- 
raphy varies as a function of space, the non-divergent form 
of the continuity equation involves the vertically integrated 
transport so that 

_ __ 1 d\p _ 1 d\j/ 

In this case the geostrophic statement is expressed by 
I = gVf 

where \[/ is now the transport stream function. Taking a curl 

f 

of this yields the condition that V — x = 0, which shows 

H f 

that the transport is parallel to the contours of — = constant. 

H 

In either case of H = constant or a variable it is clear that 
the stream lines and hence the flow will not be parallel to 
the contours of sea surface topography on a rotating earth 
with variable Coriolis parameter. Hence the hypothesis that, 
given a surface topography, one can deduce the magnitude 
and the direction of the near surface geostrophic current is 
not precisely valid over the global oceans in general and more 
particularly so in the low latitudes. 

In the context of altimetry, if the sea surface topography 
is given over a domain, equation (2.6) may be viewed as 
an inhomogeneous elliptic equation for the stream function 
\j/. Since the ocean basin is considered to be closed, the 
boundary condition on the stream function is simply 

\p — 0 on the boundary. (2.7) 

Equation (2.6) may then be solved in a straight forward 
manner for the stream function. Note that the stream func- 
tion involved here is not just the geostrophic stream func- 
tion but is associated with the full dynamics of equation (2.1). 
However, as mentioned earlier, the altimetric measurements 
over an ocean basin suffer from contamination due to vari- 
ous sources of error, the geoid being the main contributor. 
Hence it is not advisable to make an attempt to solve equa- 
tion (2.6) directly without subjecting the altimetric residu- 
als to a selective elimination of bad data points, and 
smoothing and extrapolating the remainder to define the sur- 
face topography at all grid points in the basin. This requires 
the adoption of some specific procedure. The most suitable 
one will be a method using some functions explicitly reflect- 
ing the properties of the basin under consideration. However, 
at this stage it is not obvious how to construct these functions. 

In the following a procedure is described that would first 
establish a set of appropriate characteristic functions for the 


stream function. Then from these stream function modes, 
a separate set of characteristic functions are derived that 
would be relevant to the analysis of surface topography. If 
the vorticity field r?(x, y) is assumed to be kndwn, equation 
(2.5) represents an inhomogeneous equation for \p. With the 
boundary condition (2.7), conceptually the problem can again 
be solved to obtain the stream function. In principle it is, 
of course, not possible to get a field of vorticity from meas- 
urements and particularly not from altimetric retrievals. A 
convenient way to look at the problem is simply to consider 
expanding the stream function in terms of the characteristic 
functions of the operator V 2 ; that is, in terms of the solu- 
tion to the problem: 

^ V'a fia 'Pa *2 8) 

\f/ a = 0 on the boundary. v ' } 

p a " s are the characteristic values associated with the 
vectors \[/ a . It is clear that (2.8) is a self-adjoint problem and 
the iAc/s form an orthogonal set. The scalar index a is a 
proxy to an ordering of the i/^’s in some as yet unspecified 
manner. Since the i/^’s constitute a complete set, they are 
capable of representing any arbitrary field over the basin. 
Hence we can express 

'P = XP a 'Pa ( 2 . 10 ) 

a 

Let the orthogonality of the i/^’s be defined by 

1 'Pa 'Pe dA. = , (2.11) 

where d aj3 is the Kronecker delta. The integration is per- 
formed over the area of the basin. The expansion coefficients 
are then given by 

P a = l^a \ 'P'Pa dA. (2.12) 

When the sum on the right hand side of equation (2.10) 
spans the complete spectrum of the elliptic operator in (2.8), 
equation (2.10) represents a least squares approximation to 
p with the usual restrictions on quadratric integrability and 
continuity of ^ and its derivatives. The stream function in 
(2.10) when it is substituted into the governing dynamical 
equations, determines the velocity field generated by the sur- 
face wind stress in the ocean basin under consideration. Con- 
vergence of the series (2. 10) is determined by the expansion 
coefficients P a , whose value is governed by the properties 
of the imposed wind field. The kinetic energy of the circu- 
lation in the basin is given by 

KE = V 2 ( H(u 2 + v 2 ) dA - Vt E HP„ 2 . 

J a 

Hence the contribution from each of the “modes” to the 
variance of the circulation is governed by the wind forcing 
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(2.18) 


and presumably the identity of the dominant modes would 
change with different wind fields. Consequently, in carry- 
ing out the objective analysis, it will be useful to have an 
a priori idea of the dominant modes, say for example, as 
a function of the seasonal fields. Note that the shape of the 
basin influences the properties of the stream function modes 
through equation (2.8). If real topography is also taken into 
account that would influence the shapes of the modes through 
a variant of (2.8) as considered in Rao and Schwab (1981). 

Substitute the expansion (2.10) into the vorticity equa- 
tion (2.4) and use the orthogonality condition (2.11). The 
result is: 


^- + (xn-c)P = V, 


(2.13), 


where we have defined the column vectors 7 and 7 


7 = Column (P„) 

(2.14) 

F - Column (FJ - - i j dA. 

In (2.13), II is the identity matrix and the elements the 
matrix C are defined by 

C s Q = U^ dA. (2.15) 

Hence, if the wind stress is prescribed, equation (2.13) 
can be solved directly for the expansion coefficients P a and 
the stream function and the velocity fields can be obtained 
through (2.10, 2.3). If in particular, the wind stress is steady 
and only the steady state solution is sought, it is simply given 
by 


V = (XII - C)- 1 V (2.16) 

Having obtained \ p and the flow field, it remains to es- 
tablish a basis to compute the associated height field solu- 
tion. This is accomplished through the use of the divergence 
equation (2.6) which may be viewed now as an inhomogene- 
ous equation for f . A convenient way to develop the height 
field solution is to express it in terms of the following ex- 
pansion: 


r = 5 r« (2.i7) 

f a ’ s are characteristic functions for the height field. Note 
that in the above representation for the f field the expansion 
coefficients are deliberately chosen to be the same as those 
for the stream function given in (2.10). In view of the fact 
that the expansion coefficients in (2. 10 and 17) are the same 
substituting the above in (2.6) one obtains 


gV 2 f a = V • 

which is simply a balance equation establishing a relation 
between the characteristic modes of the height field and those 
of the stream function field. The operators involved on both 
sides of the equation are elliptic. The important aspect to 
note is that the zonal dependence of the operators is the same 
on both sides of the equation but the latitudinal dependence 
differs because of the variation of the Coriolis parameter in 
that direction. In view of this the boundary conditions for 
are: 


f a = 0 on meridional boundaries 
3f a dxl/ a 

g -r— = f -r— on zonal boundaries 
& dy dy 


(2.19a) 

(2.19b) 


The second condition follows from the second momen- 
tum equation of (2. 1) using the assumption of zero wind stress 
normal to the zonal boundaries. This assumption permits the 
determination of basis functions for the height field once for 
all for any given basin just like the stream function modes. 
Since the primary wind stress force for global oceans is 
predominantly zonal, the assumption of zero meridional wind 
stress on the north-south boundaries is certainly a very 
reasonable one. Equation (2.18) subject to the conditions 
(2.19) can be solved by numerical methods for any natural 
basin. The f„’s substituted into equation (2. 17) give the sur- 
face height field associated with the wind driven circulation. 

In the above process a set of characteristic functions for 
the height field have been established that would be useful 
for objective analysis of altimetric residuals as will be shown 
later. In addition, since the expansion coefficients for both 
the height and flow fields are the same, once the P a ’s are 
calculated through objective analysis of the altimetric meas- 
urements, the stream function field associated with the dy- 
namics of equation (2. 1) is also obtained. As is evident the 
flow field associated with these dynamics is the total flow 
field and not just the geostrophic component corresponding 
to the surface topography. 

OBJECTIVE ANALYSIS 

Suppose that altimetric residuals describing the sea sur- 
face topography are available at certain discrete number of 
points over a given basin. In view of the uncertainties in- 
volved in these data due to the sources of error mentioned 
earlier, some of the data would get rejected selectively. Ac- 
ceptable data points obtained from this process would gener- 
ally not be given over a regularly spaced grid. The resulting 
sea surface height field, as shown in equation (2.17), can 
be represented spectrally by: 

r = £ p a r„ (3.D 

a - 1 
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The series is truncated at N modes and the N coefficients 
P a are the unknowns that need to be determined from the 
data. If the data are available at M number of points and each 
point is designated as 

fj = f(Xj,yj); j = l, 2. . .M, 
then equation (3.1) becomes 

fj = E P« fa(Xj,yj). (3-2) 

J a = 1 J J 


Equation (3.2) represents M number of equations in N 
unknowns. If M = N, a set of N equations in as many 
unknowns are obtained. The expansion coefficients thus de- 
termined would exactly fit the data at the observation points 
leading to a zero rms error at the data points. In general, 
however, it is preferable to take N < M so that there are 
fewer unknowns than equations (an overdetermined system) 
and determine the coefficients in a least squares sense. That 
is, since equation (3.2) can not be satisfied exactly, we try 
to minimize the rms difference which is defined by 

r- -d/2 


1 M / N \ 2 

E = - = W 


(3.3) 


where are the height field data from altimetric measure- 
ments. If a column vector of length M x 1 is defined as 

T = Column (f|, £m) 


then equation (3.2) can be expressed in a matrix form 

T = I V (3.4) 


where B is a rectangular matrix of dimensions MxN with 
elements given by 


the basin, the objectively analyzed sea surface topography 
can be mapped over the entire basin from equation (3.1). 
As mentioned earlier, the basis functions reflect the geometric 
peculiarties of an individual basin or global oceans and hence 
would constitute a more suitable set of functions for analyz- 
ing data than any analytical functions. Also since the basis 
functions used in the objective analysis of the height field 
are derived from the solution of the inhomogeneous problem 
(2.17, 18, 19), the expansion coefficients P a can now be 
used to synthesize the stream function field over the basin 
through equation (2. 10). Hence the height field analysis tech- 
nique developed here permits the possibility of computing 
not just the geostrophic flow associated with the measured 
topography but the actual (or total) flow associated with the 
topography, Coriolis forces, surface wind stress, and bot- 
tom friction. Another aspect of the technique to notice is that 
the characteristic modes of the stream function and the height 
field are completely independent of the values assigned to 
the depth H of the wind driven ocean or the friction coeffi- 
cient in the momentum equations (2.1). Consequently these 
parameters, which are by no means well established cons- 
tants, have no influence on the analysis of the surface topog- 
raphy and the synthesis of the stream function field. 

In the frame work of the simple dynamics expressed in 
equations (2.1,2) and the analysis technique developed here, 
it is possible to make inferences on the wind stress acting 
on the ocean surface. Once the expansion coefficients P a are 
determined from the objective analysis of the altimetric 
residuals, one can immediately calculate 

7 = (XII - C) 7 (3.6) 

from equation (2.16). The elements of the column vector 7 
are essentially the projections of the wind stress curl on the 
various stream function modes as shown in equation (2. 14). 
Let Q be the stream function of the wind stress so that 


dQ 


da 


Tx dy ,Ty dy 


(3.7) 


Then we can write from equation (2.14) 


1 


F “ = --g-U«v 2 nd A 


(3.8) 


If 0 is now expressed as 


B = (Bj a ) = r a (Xj, yj) 

and V is a column vector of length N X 1 given by 

V = Column (P b P 2 , . . . . Pn) 

The elements of the matrix B = (B ja ) are the values of 
the characteristic modes at the locations j and are avail- 
able from the height field characteristic functions previous- 
ly determined. The least squares solution to equation (3.4) 
is then given by 

V = (B t B)- 1 B T ? (3.5) 

where B T is the transposed matrix. Having determined ~p , 
since the functions are known over the whole domain of 


a = £ R« (3.9) 

a 

and (3.9) is substituted into equation (3.8), the result is 
R a = H F„ 


{ 
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and 


(3.10) 


Q = L H F a 

after using equation (2.8) and the orthogonality condition 
(2.11). This step determines the stream function for the wind 
stress Q through (3.10). Therefore, it is now possible, in prin- 
ciple, to compute the rotational part of the wind stress both 
in magnitude and direction using altimetric residuals of sea 
surface topography and the analysis technique developed 
here. However, it should be noted that in the calculation of 
V and R a through equations (3.6 and 10) the 
numerical values assigned to the depth H and friction coeffi- 
cient enter explicitly. Since neither of these parameters are 
in the nature of universal constants, the inference on the mag- 
nitude of the stress would be doubtful since they directly in- 
fluence R a through (3.6 and 10). But the directionality of 
the wind stress should be more reliable since X occurs only 
in the matrix in equation (3.6) that determines F a . 

Since an altimetric mission provides wind speed but not 
direction (Brown etal, 1981) the procedure advocated here 
raises the possibility of being able to assign directionality 
to the rotational part of the mean wind stress. Such a capa- 
bility gives an opportunity to check the wind stresses meas- 
ured from an altimeter and a scatterometer for consistency. 
Even though it is a reasonable assumption, it should be borne 
in mind that the present analysis procedure is based upon 
the premise that the large scale wind stress is nondivergent. 
In any case, since the wind speeds can be obtained directly 
from the altimetric measurements it is also possible to make 
an estimate of the depth H of a “wind driven ocean” using 
equations (3.7 and 10). 

RESULTS FOR AN IDEAL BASIN CASE 



m a and n a (integers) 

The above solution for \[/ a takes into account the ortho- 
normality condition (2.11). The basic stream function modes 
are simply gyres. The number of gyres inside the basin de- 
pends on the integer values of m a and n a . Figure 1 (top left 
and right panels), computed from (4.1) on a grid of Ax = 
L/50 and Ay = B/25 shows the patterns for two modes: 
m a = 1 = n a and m a = 2 = n a . For natural basins the 
important starting point, namely the solution to (2.8), will 
have to be constructed using numerical methods. This was 
done by using a grid of Ax = L/50 and Ay = B/30 in solv- 
ing for the characteristic values and functions resulting in 
1421 interior points for 0 a . The matrix of 1421 x 1421 was 
solved using a Lanczos procedure in which each Lanczos 
function at an iteration was orthogonalized with respect to 
all the preceding ones. The details of the general Lanczos 
procedure are not discussed here for the sake of brevity (see 
Platzman 1975 for details). Figure 1 (bottom panels) show 
the numerically calculated characteristic functions that cor- 
respond to the analytic ones shown in the top panels of Figure 
1 . A comparison of the analytical and numerical solutions 
shows an excellent ability of the numerical procedure to com- 
pute the necessary basis functions needed for the present anal- 
ysis. Using (4. 1) the elements of the matrix C are given by: 

C = 4 0 m 7 m o 

7 “ L yfiT a n y (m2 - m2 ) 

provided (4.2) 

m a - m 7 is odd and n a = n 7 


In this preliminary study, the above analysis will be ap- 
plied to an ocean basin of an ideal shape. The specific case 
considered is a rectangular ocean on a 0-plane. The reason 
for considering a simple case like this is simply because the 
solution to this case is well known and the ability of the above 
objective analysis procedure using the characteristic func- 
tions in reproducing the known solution can be easily tested 
before attempting to apply the technique to more complicated 
cases of real ocean basins. In addition the problem can be 
solved by semi-analytic procedure as well as by purely nu- 
merical approach allowing the possibility to judge the effects 
of discretization errors. Consider a rectangular basin 
O < X < L and O < Y < B. The solution to the charac- 
teristic value problem (2.8) is 


The wind stress acting on the surface is taken to be 

t* = ~ Sin y Cos Y 

r y = ® Cos ^ Sin 2 (4.3) 

ft = — Sin ^ Sin ^ 

7T L B 

This wind stress field is free of divergence and represents 
a simple anticy clonic gyre, fl is the stream function associated 
with the above stress. The forcing elements F a defined in 
(2.14) are given by: 


K 


2 

VBL/v 


m„ xx n a xy 
Sin Sin -2— ^ 


L B 


(4.1) 


BVBL^ 

2xH 


for m a = n a = 1 only 


(4.4) 
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In view of the simple form of the assumed wind stress 
field, it has a non-trivial projection on only the stream func- 
tion mode that corresponds to the integers m a = 1 = n a . 
The expansion coefficients P a required in equation (2.10) 
to determine the stream function associated with steady cir- 
culation are obtained by solving equation (2.16). This can 
now be accomplished using the coefficients given in equa- 
tions (4.2 and 4). The parameter values used in the numeri- 
cal computations are: 


of Figure 2. Once again the agreement between the analytic 
functions and the numerical ones is excellent. 

Substitution of the P a ’s obtained from the solution of 
equation (2. 16) into (2.10) yields the stream function of the 
circulation driven by the wind field given in equation (4.3). 
Figure 3 (bottom) shows the circulation field synthesized 
using 60 characteristic modes. In view of the conditions stated 
in (4.2, 4) the ordering of m a and n a with respect to the sca- 
lar index a was prescribed as 


L = 10 9 cm. 

B = 2 7T .10 s cm. 
i3 = 2.315 . 10~ 13 /cm. sec 
g =981 cm/sec 2 
X = 10 _5 /sec 
H = 2.10 4 cm. 


Several values of X were tried. For X < < 10 -5 , the in- 
version of the matrix in equation (2.16) does not converge 
well. For values of X > > 10 -5 , the flow field tends to be 
very damped. As a compromise, the value of X=10 -5 was 
used in the calculations reported here. As mentioned earli- 
er, the exact values of X and H are not definable but also, 
as previously noted, the basis functions for the height and 
stream fields used in the objective analysis are not affected 
by the choice of the numerical values of these two 
parameters. Consider now the problem defined by equations 
(2.18, 19) which give the characteristic functions for the 
height field. The solution for these functions is: 


fa = 


gVBL^ a 


n a 7ry 

f Sin — + 


B 


^ B « _ I W 

Cos 


(4.5) 


where 


a + b 2 


Sin 


m^xx 


L a = m a 7r/L and B a = n a 7r/B. 


Figure 2 (top panels) show two of these height field basis 
functions corresponding to m a = l=n a and m a =2=n a . 
(Here f is simply taken as /3y .) The number of gyres is again 
determined by the integer values of m a and n a . In the east- 
west direction, the number of gyres is the same as the in- 
teger value of m a . If the integer value of n a is i, there are 
i-1 closed gyres in the north-south direction and the southern 
most gyre exhibits the pattern of an open oval. When the 
numerically determined stream function modes are used to 
solve for the height field modes from (2.18, 19), the problem 
is reduced to inverting a matrix of 1519 x 1519 using the 
51x31 grid. The matrix is very sparse and can be easily 
inverted using available computer subroutines. Figure 2 (bot- 
tom panels) shows the height field functions obtained numer- 
ically that correspond to the analytical ones in the top panels 


m a = a for a = 1, 2, . . . . N 
n a = 1 for all a 


where N is the maximum number of modes taken into con- 
sideration.* The figure shows the classical westward inten- 
sification of the circulation as shown by Stommel (1948). 
The maximum value of the transport occurs approximately 
at a distance of one fifth of the length of the basin from the 
western boundary. The precise location of the maximum is 
governed by the choice of X. For example, when 
X=0.4x 10" 5 , the maximum occurs at a distance of about 
l/10th of L from the western boundary. The vorticity equa- 
tion (2.4) can be solved exactly for the forcing function given 
by (4.3). The solution in terms of the stream function is: 


* - K + l) Sin — + Cos — 

il3L\B 2 / L L 


+ t 


1 + e 


2X 


+ $ L 


2\ 


- L- *r + 0)L 


- e 


2X 


1 + e 


- £ + 0 )L 


2X 


2X 


- 1&-9 L- £ + L 


- e 


2X 




(4.6) 


2X 


+ e x 




K = B2 ± L2 
\ir + B 


1 


B ! ) rt ! / Ajr 


*An additional 50 modes with randomly assigned integer 
values to m a and n a were included in the computation to see 
their impact on the numerical process of calculating P a . The 
calculations were completely unaffected by the presence of 
these additional modes and the P a values for all of them 
were calculated to be zero. 
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e = 


B 2 4\ 2 / 


Figure 3 (top) shows a plot of the above exact solution 
to the problem. Comparison of top and bottom panels of 
Figure 3 shows an excellent agreement between the spectral 
solution and the exact solution. 

The sea surface topography associated with the circula- 
tion can now be computed using the expansion coefficients 
P a in equation (2.17) and the height field basis functions 
given in (4.5). Figure 4 illustrates this topography which also 
exhibits a highly asymmetric pattern. Basically it consists 
of a high pressure pattern over most of the basin with the 
maximum occurring approximately at a distance of L/5 from 
the western boundary and slightly greater than B/2 from the 
southern boundary of the basin. Contours of surface topog- 
raphy intersect the equator as shown in the figure. The solu- 
tion in Figure 4 can be viewed as the sea surface topography 
obtained from a perfect altimetric mission— namely, one 
without any contamination from the various sources of er- 
ror mentioned earlier. Application of simple geostrophic con- 
siderations to compute the flow from the topography in 
Figure 4 obviously is going to provide a circulation pattern 
that will not resemble the actual one shown in Figure 3 over 
most of the basin and particularly along the southern bound- 
ary where f=0. 


In performing the objective analysis, we consider the data 
field shown in Figure 4 as the “true” data field. If one sup- 
presses a certain amount of data and performs the objective 
analysis on the remainder to reconstruct the height field, the 
“true” solution would serve as a basis to estimate the ac- 
curacy or fidelity of the analysis procedure. If the basin con- 
sidered were a real ocean basin, the elimination of the 
residuals would take into account the errors expected in the 
retrievals from the sources mentioned earlier over various 
parts of the ocean basin. Since the basin under considera- 
tion here is a simple idealized one, we have suppressed the 
data on a purely random basis and performed objective anal- 
ysis on the remaining data without introducing any errors 
in the data. The height field shown in Figure 4 was calcu- 
lated on a grid of 51 x26 which gives a total number of 1274 
data points because of the boundary condition (2.19). The 
analysis was performed on randomly selected sets of 300, 
600, and 900 error free data points using in each case 30 
and 40 modes. The results are summarized in Tables 1 and 2. 

An examination of Tables 1 and 2 shows clearly that the 
objective analysis has done a remarkable job in recovering 
the “true” solution with as few as approximately one quarter 
of the total data points available using 30 modes. Figure 5 
shows the surface topography (top panel) and the stream func- 
tion field (bottom panel) recovered from the objective anal- 
ysis using 30 modes and 300 randomly selected data points. 


2-D WIND FIELD - 110 MODES (MIXED) 
FINAL ZETA PLOT 
LAMBDA =0.100E-0 



1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49 51 


Figure 4: Sea surface topography pattern (in cms.) associated with the circulation shown in 
Figure 3. 
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2-D ANALYZED ZETA 
M =30 N = 300 






Table 1 : Results of objective analysis for Z, (in cm) using error free data. 


No. of 
modes 

No. of 
points 

RMS difference 
for Zeta at data 
points 

RMS difference 
for Zeta over the 
basin 

Maximum and 
minimum value 
for Zeta from 
analysis 

theory 

30 

300 

0.0027 

0.0017 

45.54,-17.54: 

45.53, 17.54 


600 

0.0024 

0.0017 

45.53,-17.54 



900 

0.0019 

0.0017 

45.53,-17.54 


40 

300 

0.0007 

0.0004 

45.55,-17.54 



600 

0.0005 

0.0004 

45.55,-17.54 



900 

0.0004 

0.0004 

45.55,-17.54 



Table 2: Results for ip {in cm 2 /sec) 

No. of 

No. of 

RMS difference 

RMS difference 

Maximum 

* 

modes 

points 

for psi at data 

for psi over the 

value 




points 

basin 

for psi from 



(x 10 s ) 

( x 10 s ) 

analysis 

(xlO 9 ) 

theory 

30 

300 

0.345 

0.215 

0.5899 

0.5903 


600 

0.309 

0.213 

0.5903 



900 

0.236 

0.213 

0.5903 


40 

300 

0.088 

0.049 

0.5901 



600 

0.066 

0.048 

0.5902 



900 

0.053 

0.047 

0.5902 



*Minimum value for \f/ is zero on the boundary 


Comparison with Figures (4 and bottom panel of 3) shows 
the closeness of the analyzed patterns to the theoretical ones. 
Figure 6 shows the stream function Q of the wind stress field 
obtained through the objective analysis procedure (top panel) 
and the actual one given in (4.3) (bottom panel) that was used 
in the theoretical solution. Once again it is seen, at least in 
the simple case considered here, that the analysis procedure 
yields an unambiguous recovery of the rotational part of the 
driving force. As cautioned earlier, the magnitude of the wind 
stress depends on the numerical values assigned to X and H. 

The effect of assuming a X different from the one that 
has been used in generating the “true” solution on the 
retrieval of the wind stress stream function from the objec- 
tive analysis is shown in Figure 7. In the theoretical calcu- 
lation, X was taken to be 10 -5 . The X values used in Figure 
7 are 0.5 x 10~ 5 (top) and 10 -6 (bottom). It is seen from 
these figures that the stream function pattern of the wind 
stress is still reproduced fairly well over most of the basin. 
It is only in the westernmost part of the basin, where the 
circulation is strongest and hence is a region of largest fric- 
tion, do the actual and retrieved wind field directions differ. 
The magnitude, of course, has also changed. These differ- 


ent values of X, however, have absolutely no effect on the 
analysis for the surface topography and the stream function 
of the flow field presented in Figure 5 since the calculation 
of P c , from the data does not involve the use of X or H. 

In order to consider the effect of errors in the data on 
the objective analysis procedure, an error of + 10% was 
introduced randomly into the 300 data points. The result of 
analysis using 30 height field basis functions is shown in 
Figure 8 for the height field (top left), flow field (top right) 
and the wind stress stream function with X =0.5 x 10“ 5 
(bottom left) and 10 -6 (bottom right) instead of 10 -5 that 
was used in generating the theoretical data base. The result- 
ing fields were subjected to a smoothing according to: 

Cy = 0.8 Cy + 0.05 [C i+IJ + Ci-u + Cy + , + Cij-,] 

As expected the contours exhibit some wiggles instead of 
being smooth as with the error free data. The rms differ- 
ences and the maximum and minimum values resulting from 
the analysis of the data with noise are shown in Table 3 (note 
that these numbers are independent of the X value): 
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Table 3: Results from data with introduction of errors. 



RMS difference 

Maximum 


Minimum 


at data points 

whole basin 


of the 

dependent variable 


f (cm) 

0.025 

0.012 

45.88 


-17.55 

\p (cm 2 /sec) 

3.2 (10 5 ) 

1.5 (10 s ) 

5.93 (10 8 ) 


0 


The differences are still quite small and the maximum 
and minimum values of the height field and the stream func- 
tion of the flow field are determined very well. The patterns 
of these fields still exhibit a close resemblance to the theo- 
retical solutions. Once again, the direction of the wind stress 
resulting from the analysis agrees with the actual one over 
a large part of the basin with the exception of the western- 
most region as noted previously. 

Finally it should also be noted that as far as the objective 
analysis of the height field and recovery of the flow field 
are concerned, there is no implied assumption that the ocean 
dynamics are steady. The height field basis functions are 
related to those of the flow field through a purely diagnostic 
equation (2.18) because both the ocean circulation and the 
wind stress are assumed to be non-divergent. The expansions 
(2.10 and 17) for the forced stream function and the surface 
topography are still valid even if the circulation is non-steady. 
In the above examples, we have chosen to present the results 
of objective analysis for a steady-state situation only for con- 
venience. 


CONCLUSION 

Using a simple dynamical model of wind driven ocean 
circulation with linear bottom friction, an analytical basis 
was developed to objectively analyze the sea surface topog- 
raphy derived from altimetric residuals and, in the process, 
deduce the total flow instead of just the geostrophic compo- 
nent of the circulation. The primary assumptions of the the- 
ory are that both the ocean circulation and the surface wind 
stress are free of divergence. The method depends on con- 
structing a set of basis functions for the stream function of 
the flow field and a set of basis functions for the height field. 
The basis functions for the flow field are solutions of a 
homogeneous elliptic equation. The analysis shows that the 
appropriate set of basis functions for the height field are solu- 
tions of an inhomogeneous balance equation. When the func- 
tions are chosen in this particular way, the theory shows that 
in the spectral representatation of an arbitrary flow field and 
the associated height field in terms of these basis functions, 
the expansion coefficients for both fields are the same. 


Hence, when the sea surface topography is prescribed from 
altimetric measurements, after eliminating data in regions 
of serious uncertainty in the retrievals due to various sources 
of error, the height field functions can be used to objective- 
ly analyze the remaining data through least squares proce- 
dures and determine the expansion coefficients. The basis 
stream function modes then yield the total stream function 
of the flow associated with the topography. 

Even though the dynamical basis adopted here is a sim- 
ple wind driven circulation model of Stommel, it is certain- 
ly more complicated than the simple application of 
geostrophic dynamics in extracting ocean circulation infor- 
mation from altimetry. In particular, the method described 
here does not suffer from the fatal breakdown in the equatori- 
al regions like the geostrophic dynamics do. Considering that 
the tropical ocean circulation is one of the most important 
ingredients in the air-sea interactions on a climatic scale, the 
method outlined here allows the possibility of extracting the 
oceanographic signal in the tropics in a more reliable man- 
ner. Moreover, since the basis functions can be constructed 
for an ocean basin taking into account the geometric ir- 
regularities and the necessary boundary conditions, they 
would constitute a better set of functions in analyzing the 
data than any analytical functions that are oblivious of the 
basin geometry. An additional feature of the analysis is that 
it allows the possibility of assigning a direction to the sur- 
face wind stress that is responsible for maintaining the given 
surface topography. Even though this study considered only 
linear wind driven ocean dynamics, the method can be ex- 
tended to consider nonlinear aspects if necessary. 

The numerical modelling of ocean circulation has not yet 
matured to the stage of using real data to initialize the models 
and to consider problems of data assimilation like the at- 
mospheric models do; sparsity of oceanographic data being 
at least one of the reasons. With the advent of various ocean 
related satellite missions and in-situ data collection programs 
under planning in conjunction with the World Climate 
Research Program, substantial amounts of oceanographic 
data will soon be available. Some of the considerations 
presented in this study would become useful in initializing 
and balancing different data fields for iuture numerical ocean 
modeling efforts. 
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anc 




The theory was tested in this preliminary study on an ideal 
rectangular basin located north of the equator on a 0-plane. 
First a theoretical solution for the stream function of the flow 
and the resulting surface topography for a given wind stress 
were obtained. Using the theoretical field of the topography 
as a '‘true” data field, a randomly selected set of data points 
were subjected to the objective analysis using the height field 
basis modes through a least squares technique. The expan- 
sion coefficients thus determined and the corresponding basis 
functions were then used to reconstruct the height field and 
the flow field over the whole basin. Comparison with the 
“true” data fields shows that the objective analysis proce- 
dure does an exceedingly good job in recovering the “true” 
solution with as few as a quarter of the total available data 
in this ideal basin case. It was also shown that the effect of 
introducing random error of + 10% at all data points does 


not seriously affect the results. As pointed out in Section 3, 
the calculation of the direction of the wind stress depends 
on the value assumed for the linear bottom friction coeffi- 
cient. The calulations presented show that the direction of 
the wind stress is given with reasonable accuracy over most 
of the basin even when the bottom friction coefficient was 
varied over an order of magnitude. 

The technique discussed here is fully adaptable to com- 
plete numerical treatment and hence is applicable to real 
ocean basins. Calculations are in progress on the Atlantic- 
Indian Ocean system with the effect of topography also taken 
into account and will be reported shortly. 
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